Información

  • Author: Jose Carlos Molano de Oro

  • University: Pontificia Universidad Javeriana

  • Course: Linear Regression Analysis

  • Semester: 2022-3

  • Professor: Mario Gregorio Saavedra Rodriguez

  • Author Email: jose_molano@javeriana.edu.co

  • Professor Email: saavedrarmg@javeriana.edu.co

Paquete en R palmerpenguins

Descripción del Paquete

El paquete palmerpenguins contiene los datasets:

  • penguins: Size measurements for adult foraging penguins near Palmer Station, Antarctica
  • penguins_raw (penguins): Penguin size, clutch, and blood isotope data for foraging adults near Palmer Station, Antarctica

Se hara el análisis de regresión lineal en base al dataset penguins que posee información de:

  • especie (species)
  • isla (island)
  • sexo (sex)
  • año (year)

se tienen de igual manera para el grupo anterior contiene las medidas de

  • longitud de pico (bill_length_mm)
  • altura del pico (bill_depth_mm)
  • longitud de aleta (flipper_length_mm)
  • masa del cuerpo (body_mass_g)

Librerias en R

library(tidyverse)
library(palmerpenguins)
library(reactable)
library(reactablefmtr)
library (DT)
library(Hmisc)
library(GGally)
library(corrplot)
library(ggfortify)

Datos dataset penguins

reactable(penguins,rownames = TRUE)

Resumen dataset penguins

summary(penguins)
##       species          island    bill_length_mm  bill_depth_mm  
##  Adelie   :152   Biscoe   :168   Min.   :32.10   Min.   :13.10  
##  Chinstrap: 68   Dream    :124   1st Qu.:39.23   1st Qu.:15.60  
##  Gentoo   :124   Torgersen: 52   Median :44.45   Median :17.30  
##                                  Mean   :43.92   Mean   :17.15  
##                                  3rd Qu.:48.50   3rd Qu.:18.70  
##                                  Max.   :59.60   Max.   :21.50  
##                                  NA's   :2       NA's   :2      
##  flipper_length_mm  body_mass_g       sex           year     
##  Min.   :172.0     Min.   :2700   female:165   Min.   :2007  
##  1st Qu.:190.0     1st Qu.:3550   male  :168   1st Qu.:2007  
##  Median :197.0     Median :4050   NA's  : 11   Median :2008  
##  Mean   :200.9     Mean   :4202                Mean   :2008  
##  3rd Qu.:213.0     3rd Qu.:4750                3rd Qu.:2009  
##  Max.   :231.0     Max.   :6300                Max.   :2009  
##  NA's   :2         NA's   :2

Dado que en las columnas numéricas se tienen valores nulos (NA), estos valores se añadiran con el valor central (mediana) de cada atributo. De igual manera se expresara los años en variable de tipo factor, con fin de conocer cuantos pinguinos se tuvieron en cuenta en el estudio por año

penguins[,3:6]<-penguins[,3:6] %>% mutate_all(~ifelse(is.na(.x), median(.x, na.rm = TRUE), .x))
penguins$year<-as.factor(unlist(penguins$year))
summary(penguins)
##       species          island    bill_length_mm  bill_depth_mm  
##  Adelie   :152   Biscoe   :168   Min.   :32.10   Min.   :13.10  
##  Chinstrap: 68   Dream    :124   1st Qu.:39.27   1st Qu.:15.60  
##  Gentoo   :124   Torgersen: 52   Median :44.45   Median :17.30  
##                                  Mean   :43.92   Mean   :17.15  
##                                  3rd Qu.:48.50   3rd Qu.:18.70  
##                                  Max.   :59.60   Max.   :21.50  
##  flipper_length_mm  body_mass_g       sex        year    
##  Min.   :172.0     Min.   :2700   female:165   2007:110  
##  1st Qu.:190.0     1st Qu.:3550   male  :168   2008:114  
##  Median :197.0     Median :4050   NA's  : 11   2009:120  
##  Mean   :200.9     Mean   :4201                          
##  3rd Qu.:213.0     3rd Qu.:4750                          
##  Max.   :231.0     Max.   :6300

Valores Promedio en Columnas Muméricas (mean)

datatable(as.matrix(sapply(penguins[,3:6],function(x) mean(x, na.rm=TRUE))))

Las medidas promedio en toda la población de la muestra de pinguinos en (mm) y su masa en (g) son:

  • longitud de pico (bill_length_mm): 43.925
  • altura del pico (bill_depth_mm): 17.152
  • longitud de aleta (flipper_length_mm): 200.892
  • masa del cuerpo (body_mass_g): 4200.872

Matriz y Gráfico de Correlación entre Variables Numéricas

Matriz de Correlación

c=cor(penguins[,3:6])
y=as.data.frame(c)
y[y==1]<-" "
y <- mutate_all(y, function(x) as.numeric(as.character(x)))


reactable(as.data.frame.array((y)),
          defaultColDef = colDef(
            style = highlight_min_max(as.data.frame.array((y)))))

Gráfico de Correlación

corrplot(cor(penguins[,3:6],use = "complete.obs"),method="number")

corrplot(cor(penguins[,3:6],use = "complete.obs"),method="circle")

Se puede apreciar que los atributos que están más correlacionados positivamente entre si por orden son:

  • masa del cuerpo (body_mass_g) y longitud de aleta (flipper_length_mm): 0.87
  • longitud de aleta (flipper_length_mm) y altura del pico (bill_depth_mm): 0.66
  • masa del cuerpo (body_mass_g) y longitud de pico (bill_length_mm): 0.59

Para la construcción del modelo de regresión se construira un modelo entre las variables de masa del cuerpo (body_mass_g) vs. longitud de aleta (flipper_length_mm) y longitud de aleta (flipper_length_mm) y altura del pico (bill_depth_mm).

Para las variables masa del cuerpo (body_mass_g) vs. longitud de pico (bill_length_mm) se descarta la construcción del modelo dado que su estadistico de correlación esta muy cercano al 0.5, es decir que no se podria concluir bien si sus valores están directamente correlacionados o no.

Gráficos de Dispersión - Correlaciones de Pearson

Gráfico de Dispersión con el Tipo de Especie

ggscatmat(penguins, columns=3:6, color="species")

A partir del gráfico anterior se puede apreciar lo siguiente:

  • El gráfico de dispersión de las variables longitud de aleta (flipper_length_mm) vs masa del cuerpo (body_mass_g) tiene un comportamiento muy similar a una función lineal con pendiente positiva teniendo en cuenta todas las especies. A mayor peso del pinguino mayor sera su longitud de aleta.

  • El gráfico de dispersión de las variables longitud de pico (bill_length_mm) vs masa del cuerpo (body_mass_g) tiene un comportamiento muy similar a una función lineal con pendiente positiva teniendo en cuenta solo las especies de tipo Adelie (0.67) y Gentoo (0.59). A mayor peso del pinguino mayor sera su longitud de pico.

  • El gráfico de dispersión de las variables longitud de pico (bill_length_mm) vs longitud de aleta (flipper_length_mm) tiene un comportamiento muy similar a una función lineal con pendiente positiva teniendo en cuenta solo las especies de tipo Gentoo (0.7). A mayor longitud de pico mayor sera la longitud de aleta del pinguino.

Se proceden a construir algunos modelos de regresión lineal teniendo en cuenta las variables que tienen mas correlación con los gráficos anteriores.

Modelos de Regresión Lineal

1. Longitud de Aleta vs Masa del Cuerpo por Especie

ggplot(penguins, aes(y = flipper_length_mm, 
                      x = body_mass_g, 
                      )) +
     geom_point() +
  stat_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Ecuación de Regresión Lineal

\[ \textrm{flipper_length_mm}=\beta_0+\beta_1\times \text{body_mass_g}+e \]

\[ e_i\sim\mathcal{N}(0,\sigma^2) \] \[ f(x)=m\cdot x+b+e \]

Ajuste del Modelo

model1<-lm(flipper_length_mm~body_mass_g,data=penguins)
model1
## 
## Call:
## lm(formula = flipper_length_mm ~ body_mass_g, data = penguins)
## 
## Coefficients:
## (Intercept)  body_mass_g  
##   136.71105      0.01528

\[ \widehat{\textrm{flipper_length_mm}}_i=\widehat\beta_0+\widehat\beta_1\times \text{body_mass_g}_i+e_i \] \[ =1.367\times 10^2+1.528\times10^{-2}\times \text{body_mass_g}_i+e_i \] \[ e\sim\mathcal{N}(0,\sigma^2) \] \[ \widehat{\textrm{flipper_length_mm}}_i=\widehat\beta_0+\widehat\beta_1\times \text{body_mass_g}_i \] \[ =1.367\times 10^2+1.528\times10^{-2}\times \text{body_mass_g}_i \] \[ e_i=\textrm{flipper_length_mm}_i-\widehat{\textrm{flipper_length_mm}}_i \]

Resumen del Modelo

summary(model1)
## 
## Call:
## lm(formula = flipper_length_mm ~ body_mass_g, data = penguins)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.7543  -4.8130   0.9608   5.0946  16.6487 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.367e+02  1.990e+00   68.68   <2e-16 ***
## body_mass_g 1.528e-02  4.655e-04   32.82   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.894 on 342 degrees of freedom
## Multiple R-squared:  0.759,  Adjusted R-squared:  0.7583 
## F-statistic:  1077 on 1 and 342 DF,  p-value: < 2.2e-16

A partir de lo anterior se puede observar:

  • \(R^2 \text{ajustado}=0.7583\) el modelo tiene un buen ajuste con las variables seleccionadas para el estudio (cercano a 1).
  • \(p-\text{value}=2.2\times10^{-16} <0.05\) esto significa que la hipótesis nula del modelo se rechaza, es decir que con las variables seleccionadas existe un tipo de relación lineal entre ellas.

Representación del Modelo

ggplot(penguins, aes(x = body_mass_g, 
                      y = flipper_length_mm, 
                      )) +
     geom_point() +
  stat_smooth(method=lm, se = TRUE, color = "red")
## `geom_smooth()` using formula 'y ~ x'

Intervalos de Confianza para los Parámetros del Modelo

confint(model1)
##                    2.5 %      97.5 %
## (Intercept) 132.79589823 140.6261928
## body_mass_g   0.01436252   0.0161937

Condiciones de Aceptación del Modelo de Regresión Lineal

Residuos del Modelo
penguins$fitted_values<-model1$fitted.values
penguins$residuals<-model1$residuals

ggplot(penguins, aes(x = fitted_values, y = residuals)) +
    geom_point(aes(color = residuals)) +
    scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
    geom_hline(yintercept = 0) +
    geom_segment(aes(xend = fitted_values, yend = 0), alpha = 0.2) +
    theme_bw() +
    theme( legend.position = "none")

Los residuos se distribuyen de manera aleatoria entorno a la recta \(y=0\), por lo cual se acepta la linealidad.

Histograma de los Residuos
ggplot(penguins, aes(x = residuals)) +
    geom_histogram(aes(y = ..density..),bins=90,fill='white',color=4) +
    theme_light()

Verificación de Supuestos
autoplot(model1)

shapiro.test(model1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model1$residuals
## W = 0.98335, p-value = 0.0005215

Dado que en mi test anterior, mi \(p-\text{value}<0.05\), se rechaza la hipótesis nula por lo cual los residuos no presentan un comportamiento normal.

ggplot(penguins, aes(x = fitted_values, y = residuals)) +
    geom_point(aes(color = residuals)) +
    scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
    geom_segment(aes(xend = fitted_values, yend = 0), alpha = 0.2) +
    geom_smooth(se = FALSE, color = "firebrick") +
    geom_hline(yintercept = 0) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

2. Longitud de Pico vs Masa del Cuerpo por Especies Adelie y Gentoo

AdelieGentoo<-subset(penguins,species!="Chinstrap")
ggplot(AdelieGentoo, aes(y = bill_length_mm, 
                      x = body_mass_g, 
                     )) +
     geom_point() +
  stat_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Ecuación de Regresión Lineal

\[ \textrm{bill_length_mm}=\beta_0+\beta_1\times \text{body_mass_g}+e \]

\[ e\sim\mathcal{N}(0,\sigma^2) \] \[ f(x)=m\cdot x+b+e \]

Ajuste del Modelo

model2<-lm(bill_length_mm~body_mass_g,data=AdelieGentoo)
model2
## 
## Call:
## lm(formula = bill_length_mm ~ body_mass_g, data = AdelieGentoo)
## 
## Coefficients:
## (Intercept)  body_mass_g  
##   19.230360     0.005441

\[ \widehat{\textrm{bill_length_mm}}_i=\widehat\beta_0+\widehat\beta_1\times \text{body_mass_g}_i+e_i \] \[ =1.923+5.441\times10^{-3}\times\text{body_mass_g}+e_i \] \[ e_i\sim\mathcal{N}(0,\sigma^2) \] \[ \widehat{\textrm{bill_length_mm}}_i=\widehat\beta_0+\widehat\beta_1\times \text{body_mass_g}_i \] \[ =1.923+5.441\times10^{-3}\times\text{body_mass_g} \]

\[ e_i=\textrm{bill_length_mm}_i-\widehat{\textrm{bill_length_mm}}_i \]

Resumen del Modelo

summary(model2)
## 
## Call:
## lm(formula = bill_length_mm ~ body_mass_g, data = AdelieGentoo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.5720 -1.7594 -0.0742  1.6852  7.4499 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.923e+01  7.977e-01   24.11   <2e-16 ***
## body_mass_g 5.441e-03  1.815e-04   29.98   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.508 on 274 degrees of freedom
## Multiple R-squared:  0.7664, Adjusted R-squared:  0.7655 
## F-statistic: 898.9 on 1 and 274 DF,  p-value: < 2.2e-16
  • \(R^2 \text{ajustado}=0.7655\) el modelo tiene un buen ajuste con las variables seleccionadas para el estudio (cercano a 1).
  • \(p-\text{value}=2.2\times10^{-16} <0.05\) esto significa que la hipótesis nula del modelo se rechaza, es decir que con las variables seleccionadas existe un tipo de relación lineal entre ellas.

Representación del Modelo

ggplot(AdelieGentoo, aes(x = body_mass_g, 
                      y = bill_length_mm, 
                      )) +
     geom_point() +
  stat_smooth(method=lm, se = TRUE, color = "red")
## `geom_smooth()` using formula 'y ~ x'

Intervalos de Confianza para los Parámetros del Modelo

confint(model2)
##                    2.5 %       97.5 %
## (Intercept) 17.659872328 20.800848515
## body_mass_g  0.005083984  0.005798569

Condiciones de Aceptación del Modelo de Regresión Lineal

Residuos del Modelo
AdelieGentoo$fitted_values<-model2$fitted.values
AdelieGentoo$residuals<-model2$residuals

ggplot(AdelieGentoo, aes(x = fitted_values, y = residuals)) +
    geom_point(aes(color = residuals)) +
    scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
    geom_hline(yintercept = 0) +
    geom_segment(aes(xend = fitted_values, yend = 0), alpha = 0.2) +
    theme_bw() +
    theme( legend.position = "none")

Los residuos se distribuyen de manera aleatoria entorno a la recta \(y=0\), por lo cual se acepta la linealidad.

Histograma de los Residuos
ggplot(AdelieGentoo, aes(x = residuals)) +
    geom_histogram(aes(y = ..density..),bins=90,fill='white',color=4) +
    theme_light()

Verificación de Supuestos
autoplot(model2)

shapiro.test(model2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model2$residuals
## W = 0.99581, p-value = 0.6686

Dado que en mi test anterior, mi \(p-\text{value}>0.05\), se rechaza la hipótesis alternativa por lo cual los residuos presentan un comportamiento normal.

ggplot(AdelieGentoo, aes(x = fitted_values, y = residuals)) +
    geom_point(aes(color = residuals)) +
    scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
    geom_segment(aes(xend = fitted_values, yend = 0), alpha = 0.2) +
    geom_smooth(se = FALSE, color = "firebrick") +
    geom_hline(yintercept = 0) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Conclusiones de los 2 Modelos

Modelo 1

  • \(R^2 \text{ajustado}=0.7655\) el modelo tiene un buen ajuste con las variables seleccionadas para el estudio (cercano a 1).
  • \(p-\text{value}=2.2\times10^{-16} <0.05\) esto significa que la hipótesis nula del modelo se rechaza, es decir que con las variables seleccionadas existe un tipo de relación lineal entre ellas.
  • Teniendo en cuenta el test de normalidad de Shapiro-Wilk para los residuos del modelo, mi \(p-\text{value}<0.05\), se rechaza la hipótesis nula por lo cual los residuos no presentan un comportamiento normal. Quiza podría darse el caso que estandarizando las variables del modelo, sus residuos logren tener un comportamiento normal

Modelo 2

  • \(R^2 \text{ajustado}=0.7655\) el modelo tiene un buen ajuste con las variables seleccionadas para el estudio (cercano a 1).
  • \(p-\text{value}=2.2\times10^{-16} <0.05\) esto significa que la hipótesis nula del modelo se rechaza, es decir que con las variables seleccionadas existe un tipo de relación lineal entre ellas.
  • Teniendo en cuenta el test de normalidad de Shapiro-Wilk para los residuos del modelo, mi \(p-\text{value}>0.05\), se rechaza la hipótesis alternativa por lo cual los residuos presentan un comportamiento normal.